Grabbing SPINS gradients

## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.6     v dplyr   1.0.7
## v tidyr   1.1.4     v stringr 1.4.0
## v readr   2.0.1     v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## Warning: package 'ggseg' was built under R version 4.1.3
## Loading required package: prettyGraphs
## Loading required package: ExPosition

Read in the SPINS big table

## New names:
## * `` -> ...1
## Rows: 337120 Columns: 16
## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr  (5): ROI, Task, Subject, Network, Site
## dbl (11): ...1, grad1, grad2, grad3, grad4, grad5, grad6, grad7, grad8, grad...
## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.
## New names:
## * `` -> ...1
## Rows: 337120 Columns: 16
## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr  (5): ROI, Task, Subject, Network, Site
## dbl (11): ...1, grad1, grad2, grad3, grad4, grad5, grad6, grad7, grad8, grad...
## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.

read subject data

## Rows: 354 Columns: 38
## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr  (7): record_id, site, scanner, diagnostic_group, demo_sex, subject, ses...
## dbl (22): demo_age_study_entry, scog_tasit1_total, scog_tasit2_total, scog_t...
## lgl  (9): term_early_withdraw, has_RS_grads, has_EA_grads, RS_QC_exclude, EA...
## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.
##  [1] "record_id"                        "site"                            
##  [3] "scanner"                          "diagnostic_group"                
##  [5] "demo_sex"                         "demo_age_study_entry"            
##  [7] "term_early_withdraw"              "scog_tasit1_total"               
##  [9] "scog_tasit2_total"                "scog_tasit3_total"               
## [11] "np_composite_tscore"              "np_domain_tscore_process_speed"  
## [13] "np_domain_tscore_att_vigilance"   "np_domain_tscore_work_mem"       
## [15] "np_domain_tscore_verbal_learning" "np_domain_tscore_visual_learning"
## [17] "np_domain_tscore_reasoning_ps"    "np_domain_tscore_social_cog"     
## [19] "bprs_factor_pos_symp"             "bprs_factor_total"               
## [21] "bprs_factor_neg_symp"             "qls_total"                       
## [23] "sans_total_sc"                    "bsfs_total"                      
## [25] "subject"                          "session"                         
## [27] "n_scans_task-emp"                 "n_scans_task-rest"               
## [29] "fd_mean_task-emp"                 "fd_mean_task-rest"               
## [31] "has_RS_grads"                     "has_EA_grads"                    
## [33] "RS_QC_exclude"                    "EA_QC_exclude"                   
## [35] "RS_mo_exclude"                    "EA_mo_exclude"                   
## [37] "withdrawn_exclude"                "exclude_MRI"
##  [1] "scog_tasit1_total"                "scog_tasit2_total"               
##  [3] "scog_tasit3_total"                "np_composite_tscore"             
##  [5] "np_domain_tscore_process_speed"   "np_domain_tscore_att_vigilance"  
##  [7] "np_domain_tscore_work_mem"        "np_domain_tscore_verbal_learning"
##  [9] "np_domain_tscore_visual_learning" "np_domain_tscore_reasoning_ps"   
## [11] "np_domain_tscore_social_cog"      "bprs_factor_pos_symp"            
## [13] "bprs_factor_total"                "bprs_factor_neg_symp"            
## [15] "qls_total"                        "sans_total_sc"                   
## [17] "bsfs_total"

Check subject overlap

grad.sub <- spins_grads_wide$Subject[order(spins_grads_wide$Subject)]
behav.sub <- lol_spins_behav$subject[order(lol_spins_behav$subject)]

behav.sub[behav.sub %in% grad.sub == FALSE]
## [1] "sub-CMP0180" "sub-CMP0182" "sub-CMP0196" "sub-CMP0198" "sub-CMP0213"
## [6] "sub-ZHH0034"
grad.sub[grad.sub %in% behav.sub == FALSE]
##  [1] "sub-CMH0009" "sub-CMH0036" "sub-CMH0046" "sub-CMH0047" "sub-CMH0048"
##  [6] "sub-CMH0050" "sub-CMH0052" "sub-CMH0059" "sub-CMH0070" "sub-CMH0072"
## [11] "sub-CMH0097" "sub-CMH0104" "sub-CMH0114" "sub-CMH0118" "sub-CMH0126"
## [16] "sub-CMH0136" "sub-CMH0155" "sub-CMH0176" "sub-CMH0197" "sub-CMP0175"
## [21] "sub-CMP0178" "sub-CMP0183" "sub-CMP0202" "sub-MRC0001" "sub-MRC0008"
## [26] "sub-MRC0010" "sub-MRC0012" "sub-MRC0013" "sub-MRC0023" "sub-MRC0042"
## [31] "sub-MRC0044" "sub-MRC0046" "sub-MRC0047" "sub-MRC0051" "sub-MRC0058"
## [36] "sub-MRC0060" "sub-MRC0070" "sub-MRC0071" "sub-MRP0019" "sub-MRP0067"
## [41] "sub-MRP0079" "sub-MRP0093" "sub-MRP0099" "sub-MRP0105" "sub-MRP0116"
## [46] "sub-MRP0121" "sub-MRP0132" "sub-MRP0143" "sub-MRP0165" "sub-ZHH0014"
## [51] "sub-ZHH0017" "sub-ZHH0019" "sub-ZHH0020" "sub-ZHH0022" "sub-ZHH0023"
## [56] "sub-ZHH0033" "sub-ZHH0037" "sub-ZHH0040" "sub-ZHH0044" "sub-ZHH0045"
## [61] "sub-ZHH0047" "sub-ZHH0048" "sub-ZHH0052" "sub-ZHH0053" "sub-ZHH0054"
## [66] "sub-ZHH0057" "sub-ZHP0063" "sub-ZHP0065" "sub-ZHP0085" "sub-ZHP0093"
## [71] "sub-ZHP0094" "sub-ZHP0095" "sub-ZHP0097" "sub-ZHP0098" "sub-ZHP0100"
## [76] "sub-ZHP0103" "sub-ZHP0105" "sub-ZHP0119" "sub-ZHP0129" "sub-ZHP0133"
## [81] "sub-ZHP0134" "sub-ZHP0141" "sub-ZHP0142" "sub-ZHP0143" "sub-ZHP0144"
## [86] "sub-ZHP0155" "sub-ZHP0162" "sub-ZHP0168"
## kept subjects
kept.sub <- behav.sub[behav.sub %in% grad.sub] # 342
kept.sub[complete.cases(spins_behav_num[kept.sub,c(1:11, 17)]) == FALSE]
## [1] "sub-CMH0013" "sub-ZHP0081" "sub-ZHP0083" "sub-ZHP0118"
kept.sub <- kept.sub[complete.cases(spins_behav_num[kept.sub,c(1:11, 17)])] # 338

## grab the matching data

behav.dat <- spins_behav_num[kept.sub,c(1:11, 17)]
spins_grads_wide_org <- spins_grads_wide[,-1]
rownames(spins_grads_wide_org) <- spins_grads_wide$Subject
grad.dat <- spins_grads_wide_org[kept.sub,]

rm(spins_grads_wide, spins_grads_wide_org, lol_spins_behav, spins_behav_num, spins_grads_num)

grab some network colours

networks <- read_delim("../networks.txt", 
                       "\t", escape_double = FALSE, trim_ws = TRUE) %>%
  select(NETWORK, NETWORKKEY, RED, GREEN, BLUE, ALPHA) %>%
  distinct() %>%
  add_row(NETWORK = "Subcortical", NETWORKKEY = 13, RED = 0, GREEN=0, BLUE=0, ALPHA=255) %>%
  mutate(hex = rgb(RED, GREEN, BLUE, maxColorValue = 255)) %>%
  arrange(NETWORKKEY)
## Rows: 718 Columns: 12
## -- Column specification --------------------------------------------------------
## Delimiter: "\t"
## chr (4): LABEL, HEMISPHERE, NETWORK, GLASSERLABELNAME
## dbl (8): INDEX, KEYVALUE, RED, GREEN, BLUE, ALPHA, NETWORKKEY, NETWORKSORTED...
## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.
networks
## # A tibble: 13 x 7
##    NETWORK              NETWORKKEY   RED GREEN  BLUE ALPHA hex    
##    <chr>                     <dbl> <dbl> <dbl> <dbl> <dbl> <chr>  
##  1 Visual1                       1     0     0   255   255 #0000FF
##  2 Visual2                       2   100     0   255   255 #6400FF
##  3 Somatomotor                   3     0   255   255   255 #00FFFF
##  4 Cingulo-Opercular             4   153     0   153   255 #990099
##  5 Dorsal-Attention              5     0   255     0   255 #00FF00
##  6 Language                      6     0   155   155   255 #009B9B
##  7 Frontoparietal                7   255   255     0   255 #FFFF00
##  8 Auditory                      8   250    62   251   255 #FA3EFB
##  9 Default                       9   255     0     0   255 #FF0000
## 10 Posterior-Multimodal         10   177    89    40   255 #B15928
## 11 Ventral-Multimodal           11   255   157     0   255 #FF9D00
## 12 Orbito-Affective             12    65   125     0   168 #417D00
## 13 Subcortical                  13     0     0     0   255 #000000

get row and column designs

## match ROIs to networks
ROI.network.match <- cbind(spins_grads$ROI, spins_grads$Network) %>% unique
ROI.idx <- ROI.network.match[,2]
names(ROI.idx) <- ROI.network.match[,1]
### match networks with colors
net.col.idx <- networks$hex
names(net.col.idx) <- networks$NETWORK

## design matrix for subjects
sub.dx <- spins_dx_org[kept.sub,]

diagnostic.col <- sub.dx$diagnostic_group %>% as.matrix %>% makeNominalData() %>% createColorVectorsByDesign()
rownames(diagnostic.col$gc) <- sub(".","", rownames(diagnostic.col$gc))

## design matrix for columns - behavioral 
behav.dx <- matrix(nrow = ncol(behav.dat), ncol = 1, dimnames = list(colnames(behav.dat), "type")) %>% as.data.frame

behav.col <- c("scog" = "#F28E2B",
               "np" = "#59A14F",
               "bsfs" = "#D37295")

behav.dx$type <- sub("(^[^_]+).*", "\\1", colnames(behav.dat))
behav.dx$type.col <- recode(behav.dx$type, !!!behav.col)

## design matrix for columns - gradient
grad.dx <- matrix(nrow = ncol(grad.dat), ncol = 4, dimnames = list(colnames(grad.dat), c("gradient", "ROI", "network", "network.col"))) %>% as.data.frame

grad.dx$gradient <- sub("(^[^.]+).*", "\\1", colnames(grad.dat))
grad.dx$ROI <- sub("^[^.]+.(*)", "\\1", colnames(grad.dat))
grad.dx$network <- recode(grad.dx$ROI, !!!ROI.idx)
grad.dx$network.col <- recode(grad.dx$network, !!!net.col.idx)

## get different alpha for gradients
grad.col.idx <- c("grad1" = "grey30",
                    "grad2" = "grey60",
                    "grad3" = "grey90")
grad.dx$gradient.col <- recode(grad.dx$gradient, !!!grad.col.idx)

Run PLS-C

pls.res <- tepPLS(behav.dat, grad.dat, DESIGN = sub.dx$diagnostic_group, make_design_nominal = TRUE, graphs = FALSE)

## swith direction for dimension 3
pls.res$TExPosition.Data$fi[,3] <- pls.res$TExPosition.Data$fi[,3]*-1
pls.res$TExPosition.Data$fj[,3] <- pls.res$TExPosition.Data$fj[,3]*-1
pls.res$TExPosition.Data$pdq$p[,3] <- pls.res$TExPosition.Data$pdq$p[,3]*-1
pls.res$TExPosition.Data$pdq$q[,3] <- pls.res$TExPosition.Data$pdq$q[,3]*-1
pls.res$TExPosition.Data$lx[,3] <- pls.res$TExPosition.Data$lx[,3]*-1
pls.res$TExPosition.Data$ly[,3] <- pls.res$TExPosition.Data$ly[,3]*-1
PlotScree(pls.res$TExPosition.Data$eigs)

Results

Dimension 1

lxly.out[[1]]

gridExtra::grid.arrange(bar.grad1, bar.grad2, bar.grad3, ncol = 1)

PrettyBarPlot2(pls.res$TExPosition.Data$pdq$p[,1],
               threshold = 0, 
               color4bar = behav.dx$type.col, 
               horizontal = FALSE, main = "Loadings - behavioural")
## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.

Dimension 2

lxly.out[[2]]

gridExtra::grid.arrange(bar.grad1, bar.grad2, bar.grad3, ncol = 1)

PrettyBarPlot2(pls.res$TExPosition.Data$pdq$p[,2],
               threshold = 0, color4bar = behav.dx$type.col, horizontal = FALSE, main = "Loadings - behavioural")

Dimension 3

lxly.out[[3]]

gridExtra::grid.arrange(bar.grad1, bar.grad2, bar.grad3, ncol = 1)

PrettyBarPlot2(pls.res$TExPosition.Data$pdq$p[,3],
               threshold = 0, color4bar = behav.dx$type.col, horizontal = FALSE, main = "Loadings - behavioural")

Back into the brain

Dimension 1

## merging atlas and data by 'label'
## merging atlas and data by 'label'

Dimension 2

## merging atlas and data by 'label'
## merging atlas and data by 'label'

Dimension 3

## merging atlas and data by 'label'
## merging atlas and data by 'label'

Fancy figures

3D plot of the gradients

library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
data4plotly <- data.frame(spins_mean, grad4grp[rownames(spins_mean),])
color4plotly <- net.col.idx

fig <- plot_ly(data4plotly, x = ~grad1, y = ~grad2, z = ~grad3, color = ~Network, colors = color4plotly, alpha = 0.5, size = 5)  
fig <- fig %>% add_markers()
fig <- fig %>% layout(scene = list(xaxis = list(title = "Gradient 1"),
                      yaxis = list(title = "Gradient 2"),
                      zaxis = list(title = "Gradient 3")))
fig

Dimension 1

Dimension 2

Dimension 3